home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_2.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  4.7 KB  |  198 lines  |  [MATS/MATL]

  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.2 (Polynomial Calculus).
  9. % Section    4.2,    Introduction to Interpolation, Page 212
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % Proceed with an investigation of synthetic division for
  13.  
  14. % evaluating  Pn(x),  P'n(x)  and  ∫Pn(x)dx  where
  15.  
  16. % Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)
  17.  
  18. % It is assumed that the constant K of integration is K = 0.
  19.  
  20. % Pn(x), P'n(x) and ∫Pn(x)dx are evaluated by invoking the Matlab
  21.  
  22. % subroutine polyval(C,X) polyval(Cd,X) polyval(Ci,X), respectively.
  23.  
  24. pause % Press any key to continue.
  25.  
  26. clc;
  27. %     Several Taylor polynomial coefficient lists are 
  28.  
  29. % stored in M-files named;   zcos  zsin  ztan  zexp
  30.  
  31. % zacos  zasin  zatan  zcosh   zsinh   zsqrt   zlog
  32.  
  33. % zsqrt4   zinv   zemx2d2   zcosde   zsinde   zlogq
  34.  
  35. %     The user must determine the appropriate formula
  36.  
  37. % for  f'(x)  in order to graph  f'(x)  and  P'n(x).
  38.  
  39. % The user must also determine the formula for ∫f(x)dx
  40.  
  41. % in order to graph  ∫f(x)dx  and  ∫Pn(x)dx.
  42.  
  43. pause    % Press any key continue.
  44.  
  45. clc;
  46. % The algorithms are illustrated with the coefficients of the
  47.  
  48. % Taylor polynomial of degree  n = 11 for  f(x) = ln(1+x).
  49.  
  50. % Issue the command  zlog  to load the coefficients
  51.  
  52. % into the array  C.  The function name is loaded
  53.  
  54. % into the variable fun, the maximum degree is loaded into N.
  55.  
  56. % For graphing over [a,b] the endpoints are in  a  and  b.
  57.  
  58. % Load the Taylor coefficients.
  59.  
  60. [fun,dfun,ifun,x0,m,C,Ax] = zlog;
  61.  
  62.  
  63. pause    % Press any key continue.
  64.  
  65. clc;
  66. a = Ax(1,1);       % You can change the endpoint a.
  67. b = Ax(1,2);       % You can change the endpoint b.
  68. c = Ax(1,3);
  69. d = Ax(1,4);
  70. n = 11;            % You can change the degree n.
  71. C = flipud(C);
  72. C = C(1:n+1);
  73. C = flipud(C);
  74.  
  75. pause    % Press any key to graph f(x) and P(x).
  76.  
  77. clc; clg;
  78. h = (b-a)/200;
  79. X = a:h:b;
  80. x = X;
  81. Y = eval(fun);
  82. Z = polyval(C,X);
  83. axis([a b c d]);...
  84. plot(X,Y,'-g',X,Z,'-r');...
  85. hold on;...
  86. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  87. xlabel('x');...
  88. ylabel('y');...
  89. Mx1 = [fun,'  and  P'];...
  90. Mx2 = [Mx1,num2str(n),'(x)'];...
  91. title(Mx2);...
  92. grid;...
  93. axis;...
  94. hold off;...
  95. shg; pause    % Press any key to continue.
  96.  
  97. Mx1 = 'The function is  f(x) = ';
  98. Mx2 = 'The interval is  ';
  99. Mx3 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  100. Mx4 = 'The degree is  n = ';
  101. Mx5 = ', and the coefficients list  c  is:';
  102. clc,echo off,diary output,...
  103. disp(''),disp([Mx1,fun]),...
  104. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  105. disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
  106. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])'); end,...
  107. diary off, echo on
  108.  
  109.  
  110.  
  111.  
  112.  
  113. pause    % Press any key to graph f`(x) and P`n(x).
  114.  
  115. clc; clg;
  116. a = Ax(2,1);
  117. b = Ax(2,2);
  118. c = Ax(2,3);
  119. d = Ax(2,4);
  120. x = X;
  121. Y = eval(dfun);
  122. Dv = zeros(n,1);
  123. for j=1:n,
  124.   Dv(j) = n-j+1;            % Form vector [n,n-1,...,1]
  125. end
  126. Cd = C(1:n).*Dv;            % Coefficients for Pn`(x).
  127. Z = polyval(Cd,X);
  128. axis([a b c d]);...
  129. plot(X,Y,'-g',X,Z,'-r');...
  130. hold on;...
  131. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  132. xlabel('x');...
  133. ylabel('y');...
  134. Mx1 = [dfun,'  and  P`'];...
  135. Mx2 = [Mx1,num2str(n),'(x)'];...
  136. title(Mx2);...
  137. grid;...
  138. axis;...
  139. hold off;...
  140. shg; pause    % Press any key to continue.
  141.  
  142. Mx1 = 'The function is  f`(x) = ';
  143. Mx2 = 'The interval is  ';
  144. Mx3 = 'P`n(x) = d(1)x^n + d(2)x^(n-1) + ... + d(n)x + d(n+1)';
  145. Mx4 = 'The degree is  n = ';
  146. Mx5 = ', and the coefficients list  Cd  is:';
  147. clc,echo off,diary output,...
  148. disp(''),disp([Mx1,dfun]),...
  149. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  150. disp(Mx3),disp([Mx4,num2str(n-1),Mx5]),...
  151. for i=1:5:n, disp(Cd([i:min(i+4,n)])'); end,
  152. diary off, echo on
  153.  
  154.  
  155.  
  156.  
  157. pause    % Press any key to graph ∫f(x)dx and ∫Pn(x)dx.
  158.  
  159. clc; clg;
  160. a = Ax(3,1);
  161. b = Ax(3,2);
  162. c = Ax(3,3);
  163. d = Ax(3,4);
  164. x = X;
  165. Y = eval(ifun);
  166. Iv = zeros(n+1,1);
  167. for j=1:n+1,
  168.   Iv(j) = 1/(n-j+2);            % Form vector [1/(n+1),...,1/2,1]
  169. end
  170. Ci = C(1:n+1).*Iv;              % Coefficients for ∫Pn(x)dx.
  171. Ci = [Ci;0];
  172. Z = polyval(Ci,X);
  173. axis(Ax(3,:));...
  174. plot(X,Y,'-g',X,Z,'-r');...
  175. hold on;...
  176. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  177. xlabel('x');...
  178. ylabel('y');...
  179. Mx1 = [ifun,'  and  ∫P'];...
  180. Mx2 = [Mx1,num2str(n),'(x)dx'];...
  181. title(Mx2);...
  182. grid;...
  183. axis;...
  184. hold off;...
  185. shg; pause    % Press any key to continue.
  186.  
  187. Mx1 = 'The function is  ∫f(x)dx = ';
  188. Mx2 = 'The interval is  ';
  189. Mx3 = '∫Pn(x)dx = a(1)x^n + a(2)x^(n-1) + ... + a(n)x + 0';
  190. Mx4 = 'The degree is  n = ';
  191. Mx5 = ', and the coefficients list  Ci  is:';
  192. clc,echo off,diary output,...
  193. disp(''),disp([Mx1,ifun]),...
  194. disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
  195. disp(Mx3),disp([Mx4,num2str(n+1),Mx5]),...
  196. for i=1:5:n+2, disp(Ci([i:min(i+4,n+2)])'); end,
  197. diary off, echo on
  198.